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^H ■ Abstract 

^fi • No Runge-Kutta method can be energy preserving for all Hamiltonian systems. But 

for problems in which the Hamiltonian is a polynomial, the Averaged Vector Field (AVF) 
- method can be interpreted as a Runge-Kutta method whose weights bi and abscissae d 

.^^ ' represent a quadrature rule of degree at least that of the Hamiltonian. We prove that when 

►-_ \ the number of stages is minimal, the Runge-Kutta scheme must in fact be identical to the 

^ . AVF scheme. 

^: 

2 ■ 1 Introduction and main result 

We shall be concerned with canonical Hamiltonian systems 
^>; y' = ./-V//(y) = /(y), J=^li\\ (1) 

vn ; 

0\J , The numerical solution of problems of the this type has been treated extensively in the literature, 

CO ■ we refer to the monographs [SI E] and the references therein for details. Two of the most 

(T^ \ important properties of the system ([1]) are that the flow is a symplectic map and that the 

^D ' Hamiltonian ll{y) is preserved along any solution y(i). The circumstances under which various 

^N I numerical integrators inherit these two properties are by now fairly well understood. The focus 

. . ' in the present paper is the preservation of the Hamiltonian itself, we study integrators generating 

^ \ a sequence of approximations {y„} to the solution of ([l} such that H{yn) = H{yo) for all n > 1. 

In particular we consider what can be achieved when the Hamiltonian is polynomial and the 

«2j ', integrator is a Runge-Kutta method. For linear Hamiltonians, the resulting ODE is constant 

C^ ■ and any consistent Runge-Kutta scheme will reproduce the exact solution. If the Hamiltonian 

is quadratic, then the resulting ODE is linear, and the condition for preserving energy is that 

the stability function of the method satisfies R{z)R{~z) — 1. For polynomials of higher order it 

is not known to which extent Runge-Kutta methods can preserve the Hamiltonian. However, it 

was noted in [T3] that the Averaged Vector Field (AVF) method, defined as 

yn+l=yn + h f{{l~C)yn+^yn+l)d(, (2) 
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preserves the Hamiltonian for all problems of the form ([T|) . The AVF method has second order 
convergence. In particular, when the Hamiltonian is a polynomial, the integral can be exactly 
resolved a priori, the same result is obtained if the integral in ^ is replaced by a quadrature 
rule of sufficiently high order. This was observed in [3]. In fact, a standard linear quadrature 
formula with abscissae c — (ci, . . . , Cg)'^ and weights b = (&i, . . . , bs)'^ , results in a Runge-Kutta 
method in which the Butcher matrix is given as A — dP^ . This immediately shows that for any 
polynomial Hamiltonian system, there exist Runge-Kutta methods which exactly preserve the 
energy. Note also that any choice of quadrature rule of sufBciently high order yields the same 
approximation, the AVF method is reproduced exactly. 

As pointed out in |3] any energy-preserving integrator for ([T]) must obey all quadrature condi- 
tions, but for polynomial systems this can be relaxed. Letting the Hamiltonian be a polynomial of 
degree tti, a necessary condition for the energy to be preserved is that the quadrature conditions 
hold up to order m, or in terms of Runge-Kutta coefficients 

^6.cf"i = i, fc=l,...,m. (3) 

i 

Thus, in considering energy preserving Runge-Kutta methods for polynomial Hamiltonians of 
degree < m one may immediately restrict the focus to schemes whose coefficients satisfy ([3]). If 
771 = 2s then the smallest possible number of stages in the scheme is s the resulting abscissae and 
weights are those of the Gauss-Legendre quadrature rule. If to = 2s — 1 then the smallest possible 
number of stages is still s, but the quadrature rule is not uniquely given, although applying the 
corresponding Runge-Kutta method with A = dr^ yields the same result for all (c^, bi) satisfying 
([3]) for all fc < 2s — 1. In such a situation, we are interested in answering the question of whether 
the Butcher matrix A = dF is unique. We shall restrict our search to Butcher matrices satisfying 
the usual condition 

s 

'^aij =Ci, i = l,...,s. (4) 

We shall prove the following theorem 

Theorem 1. Let m > 3. Among all Runge-Kutta methods which exactly preserve all polynomial 
Hamiltonians of degree at most m, those with the minimal number of stages coincide with the 
AVF method 1^ when applied to such problems. The number of stages in these methods is 

L(m + 1)/2J. 

In general, there are energy preserving Runge-Kutta methods for polynomial Hamiltonian 
systems which do not coincide with the AVF-integrator. There also exist such methods of 
arbitrarily high order. Examples are easily obtained as composition methods based on the AVF- 
integrator, or by the collocation methods proposed in 10 and 7 . In the rest of this paper, we 
prove Theorem [TJ The technique we use can be summarized as follows 

1. The first step is to consider a set of conditions for energy preservation which are linear in 
the Butcher matrix A, these conditions are called the double bush conditions and may be 
thought of as a linear system M{A) — w where M is a linear map from R"^'* into R^ for 
some N to be specified. 

2. Then the rank of M is determined, the results are different in the even (m — 2s) and odd 
(to = 2s — 1) cases. A particular basis for the kernel of M is identified in each of the cases. 

3. The discretized AVF method, A — cb^ , represents a known solution and any other solution 
must be of the form A — cb'^ + N, where N is in the kernel of M. We show, by using 
certain nonlinear energy preserving conditions that such solutions require A^ = 0. 



The rest of the paper is organized as follows: In Section [5] we review some tools needed from 
the literature, and we provide general conditions for energy preservation of B-series methods 
needed in the proof. In sections three and four, we prove Theorem [1] for the cases of even and 
odd polynomial degree of the Hamiltonian respectively. 

2 Energy preservation for B-series methods 

In order to make the paper self contained, we begin this section by reviewing some general tools 
from [5] , see also [6] and [4] . General conditions for energy preservation are derived for integrators 
which possess a B-series expansion. These include the Runge-Kutta methods as a subclass. A 
good account of order theory and B-series can be found in the monographs [H [HI IH], but for 
completeness we recount briefly the main ingredients we need. 

Let T be the set of rooted trees, and for t € T we write |t| for its number of vertices. A forest 
is an unordered finite collection of trees from T, ^1^2- • -tq were each tree can appear several times, 
one may then write r = t^^ • • ■ tp'' for distinct members ti, . . . , tp, indicating that ti appears r^ 
times. The set of all forests is denoted T, and the order \t\ of a forest is the sum of the orders of 
each of its elements. An element of T is either the one-node tree •, or consists of a root to which 
a forest is attached, we use the notation t = [ti, . . . ,tq] or sometimes t = [t^Hj . . . t/] so that 
each distinct tree ti occurs r, times as a subtree of t. Sometimes we shall write B^{t) to denote 
the forest consisting of the subtrees of t. The symmetry coefficient is defined as 

ct(.) = 1, ff(i)=ri!---rp!a(ii)---a(tp). 

We recall that a large class of integrators, including in particular the Runge-Kutta methods, can 
be formally expanded into an infinite series in terms of derivatives of the vector field /, indexed 
by the set of rooted trees T . Writing yi = iphiv) where iph is the numerical flow map, we have 

yi = B{a, y)=y + ha{,)f{y) + J^a(I)/'(/)(y) + • • ■ + J^a(t)F(i)(y) + • • ■ (5) 

cr(I) a(t) 

Here a : T — > R is a method dependent coefficient map, and F{t) is the elementary differential 
corresponding to t, so that e.g. F(.) = /, F{1) = /'(/)• The B-series ^ can in many cases 
be conveniently extended to allow for pullback expansions of functions along the B-series map 
B(a,y), see for instance [5] or for non-commutative structures we refer to [T1I12). It requires the 
extension of a{t) to forests, setting for any forest t = t^ . . . t^ , ai^r) — a(ii)'"^ • • • a(i„i)'''" . One 
has for any real valued smooth function G, 



GiB{a, •)) - G + ha{,) G' f + hM'f G"{f, /) 



G + E T7ZT« WG^^^ (Fin), . . . , F{rq)), (6) 
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where we sum over forests t — T1T2 . . .Tq. If we apply (jH]) to the special case where G ^ H and 
introduce the elementary Hamiltonian, defined for every t = [ti . . . tq] as 

iJ(i)=ff(')(F(ii ),..., F(i,)), 

we conclude that 

H{B{a, y))^J2l^(f[ a(i,) ) H{t){y). (7) 

teT ^ ' \k=l ) 



The term indexed by the one-vertex tree • is interpreted as -ff^°''(F(0))(y) :— H{y). Thus to 
formally impose that H is conserved, i.e. H{y) = H{B{a,y)) for a map with B-series B{a,y) 
amounts to requiring that the the right hand side of Q minus the first term (with t = •) sums 
to zero. It is, however not so that the set of functions {H{t),t G T} is linearly independent. It 
is well-known that the dependency can be described by means of the Butcher product, defined 
between two trees u = [ui- ■ -Uq] G T and w G T as 

uo V = [uiU2 ■ ■ . UqV] G T. 

This product is non-commutative, and satisfies |wou| = |u|-|-|ti|. For elementary Hamiltonians 
one has 

H{uov)^~H{vou) (8) 

for any pair of trees u,v G T. The two trees u o v and v o u are topologically identical, v o u 
is obtained from u o v hy shifting the root one position. Conversely, any two trees ti and ^2 
which differ only by such a shift of the root can be represented as ti = u o v and t2 — v o u 
for a certain choice of u and v. This shifting of roots induces an equivalence relation on the 
set of trees by defining two trees to be equivalent if and only if one can be obtained from the 
other by zero or more root shifts. All trees in the same equivalence class clearly have the same 
number of vertices, and each equivalence class is called a free tree. The set of all free trees is 
denoted FT and those with precisely n vertices we call FT"" . The canonical projection is denoted 
■K : T ^f FT. For two trees u and v in the same equivalence class, we define k{u, v) to be the 
number of root shifts necessary to obtain v from u and k(u, u) — 0. A special role is played by 
those trees which have a factorization u o u, then ([S]) implies H{u o u) — 0. Any free tree which 
contains a member with such a factorization is called superfluous and we note that superfluous 
trees have an even number of vertices. The set of nonsuperfluous free trees will hereafter be 
denoted FT^, and FT". It follows from ([5]) that for two trees u and v in the same equivalence 
class, H{u) = (— l)''("^'")iJ(u) so, disposing of the superfluous trees for which H{t) ~ 0, we may 
rewrite ([7]) as follows 

H{B{a,y))-H{y)^Y. E /*'""'^W E L) H ^(^'-) (9) 

where t is some designated element in the equivalence class t, and where each u G 7r^^(t) is 
composed of subtrees as m = [Mi...Wg] {q depends on u). It is known that the elementary 
Hamiltonians corresponding to the set of nonsuperfluous free trees are linearly independent, and 
that leads us to the condition for energy preservation derived by Chartier et al. [5], saying 
that the innermost sum must vanish for every free tree. In fact, since the power of h in the 
above expression is |i| — 1 we need to consider trees in PT^'-+^ to obtain conditions for energy 
preservation to order n 

Theorem 2. J^ A map with B-series B{a, y) preserves energy up to order n if and only if 

E ^ ! ^ a{B^{u)) = 0, VtG U F^^ (10) 

ue7r-i(f) ^ ' k<7i+l 

Here t is a designated member of the equivalence class t. 

One may remark that all the conditions of the theorem must be satisfied in order for the 
corresponding method to be energy preserving for every Hamiltonian function H. However, in 



this article we are interested in the subclass of Hamiltonians which are multivariate polynomials 
of some prescribed degree m. For such H{y) one realizes that H{t) = whenever t contains a 
vertex with more than m emanating branches. On the other hand, one may verify that all H{t) 
corresponding to nonsuperfluous free trees with at most m emanating branches from any vertex 
form a linearly independent set when considered uniformly over all Hamiltonians of degree at 
most m. The following result is inspired by [3]. 

Theorem 3. Any consistent B-series method with coefficients a{t) which is energy preserving 
for all polynomial Hamiltonians of degree m satisfies the quadrature conditions of order k for 
1 < k < m, i.e. 

a([.''"^]) = T, l<k<m. 
k 

If the method is a Runge-Kutta method with abscissae Ci and weights bi, i = 1, . . . , s a necessary 

condition for energy preservation is 

i^^rcf^l, fc = l,...,m. (11) 

Proof. For 1 < fc < to we consider ([TU| for the free tree with fc + 1 vertices containing the bushy 
tree t — [•'^] (i.e. the tree consisting of k copies of the one- node tree as subtrees). There is only 
one other tree in the equivalence class, namely t' — [[•'°~^]]. Now a{B^{t)) = a{J') = a{»f' = 1 
for any consistent method. On the other hand a{B^{t')) — a{[J\^^^) and together with ([TU]) 
and the fact that n{t,t') = 1, a{t) = fc! and cr(i') — [k ~ 1)! we get the desired result. For 
Runge-Kutta methods it is well-known that the a{[J'^^]) is the left hand side of (|TT|) . D 

2.1 The double bush conditions 

A certain subset of the nonsuperfluous free trees will play a particular role here, these are the 
trees which yield linear conditions on the matrix A. We consider the double bush free trees 
that we denote tp^q for integers p and q in {1, 2, . . . ,m — 1}. Clearly tp^p is superfluous, and by 
symmetry, tp^q — tq^p, so one will typically require 1 <p<q<m~l. 
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Figure 1: The double bush free tree tp^q having p leaves on one side and q on the other 

For g = TO — 1 the maximal number of branches from a vertex is m. We state the resulting 
conditions for energy preservation in the following lemma. 

Lemma 4. Let {A, 6, c) be a Runge-Kutta scheme whose abscissae and weights satisfy the quadra- 
ture conditions (|11[) for 1 < k < m. Then the conditions for energy preservation imposed on the 
method by the double bush free trees tp^q, henceforth called the double bush conditions are 

pb^CP'^Ac'^ -qb'^C'^-^AcP = — —, l<p<q<m-l. (12) 

q+l p+l 

Here c = (ci, . . . , Cs) , C = diag{ci, . . . , Cg). Powers of c are defined componentwise. A partic- 
ular solution to the double bush conditions is given as 

Aavf = Cb'^ (13) 



Proof. The free tree tp q is the equivalence class containing the four distinct rooted trees: ti = 
[['"-'['"Wl t2 = [A'%h = {AA\ and u = [[.'"M.^]]]- 

'T '^ ^^ -^ 

Figure 2: The four trees in the double bush equivalence class ^2,3 
One has then K{ti,tj) — \i — i\. By the quadrature conditions, 

a(i3_(i2)) = a(.ra([.']) = ^, «(S-fe)) = «(.)^a([.n) = ^- 
whereas 

We also compute a{t\) — {p — l)!q! cr(i2) = flia) = p!?!, and cr(t4) = p!(g — 1)!. Substituting all 
this into the conditions ([TU|) we get the stated conditions. 

We finally substitute (J13p for yl in ((T2)) and use the quadrature conditions to get 

p g 1 1 



pb^CP-'cb'^c'^ - qb^C^-^cb^cP = 



{p+l){q+l) {q + l){p+l) q + 1 p+1 



u 



For a given s, the double bush conditions P^ define a linear operator M : R''^" -^ R2 (™-i)(™-2) 
acting on the set of Butcher matrices A. Generally, this operator depends on the quadrature 
coefficients(6i, Ci), as well as on m and s. However, in our case we shall always be concerned with 
quadrature formulas of the highest possible order, so that we have either m = 2s or m = 2s — 1. 
We can then also represent the abscissae and weights of the quadrature formula by means of 
a single real parameter C as follows: We assume that (ci, . . . ,Cs) are the distinct zeros of the 
polynomial Ps{x) — CPs~iix) where Pq is the qth degree Legendre polynomial relative to the 
interval [0, 1], and the weights are determined by solving ([3]) for k < s. We would thus have 
M = M{(^,s,m), but we shall restrict our attention the particular even and odd cases for m: 
M(0, s, 2s) and M{(, s, 2s — 1) respectively. For ease of notation we still denote the linear operator 
simply by M when it will be clear from the context whether we are considering the even or odd 
case. The following lemma is included without proof for future reference 

Lemma 5. In both the even and odd cases, the matrix A^i = (1 — c)b^ is in the kernel of M , 
i.e. M{{1 - c)b^) = 0, where 1 = (1, . . . , 1)^ G R^ 

It is useful to combine the double bush conditions into conditions involving arbitrary polyno- 
mials in C and c rather than the monomials used in the previous lemma, we shall let lip denote 
the linear space of polynomials of degree at most p. 

Lemma 6. Let (A, 6, c) he a Runge-Kutta scheme whose abscissae and weights satisfy the quadra- 
ture conditions (jlip for 1 < k < m. As.sume that it also satisfies the double bush conditions (J12l) . 
Let P ellp andQ e H, such that P(0) = Q(0) = 0. Then 

b'^P'{C)AQ{c) ~ b^Q'{C)AP{c) = P{1) f Q{t) dt - g(l) / P{t) dt. (14) 

Jo Jo 



Proof. Write P{z) — ^ , ap'Z^ and Q{z) — ^ , /3q'z'^ . Then, using Lemma 2] 

b^P'{C)AQ{c) - b^Q'(C)AP(c) = ^ ap./3g- (v'h^&-^Ac'^ - q' b^ C^' ~^ Ac^' 

p',q' 



n 



2.2 Some nonlinear conditions 



We shall introduce some conditions for energy preservation which are nonlinear in the Butcher 
matrix A, but that will be used in the final stage of the proof to eliminate the presence of 
elements from the kernel of M in A for an energy preserving integrator. The triple hush trees 



•^1—4^—1^ 



Figure 3: The triple bush free tree tp,r,g having p > 1 leaves on the left side, g > 1 leaves on the 
right side, and r > leaves in the middle 

yield conditions for energy preservation which are quadratic in A^ an example of such a tree is 
shown in Figure [31 Applying Theorem [51 and inserting the appropiate B-series coefficients for 
Runge-Kutta methods [8, section III. 1.1], we find the triple hush conditions 

= pb^ CP-^ AC Ac" - b^CAc" + i - rb^C'-^iAc" o Ac^) 

- b^CAcP + qb^ C"'^ AC'' Ac" (15) 

Using the same approach as in Lemma [6l we derive the following alternative version 

= b^P\C)AR{C)AQ{c) + h^Q\C)AR{C)AP{c) - P{l)b^ R{C)AQ{c) (16) 

-Q{l)b^R{C)AP{c)-h'^R'{C){AQ{c)QAP{c)) + R{l) [ P{t)dt f Q(i)dt 

Jo Jo 

where P G lip, Q e Uq, R e Ur, P(0) = (5(0) = 0, and where p < m - l,q < m - l,r < m - 2. 
The symbol signifies component wise product between two vectors. 

Finally, we include a free tree and its corresponding condition used to investigate an excep- 
tional case in Section [4l see Figure [4l The energy preservation condition corresponding to this 
tree is found to be 

b^iAc)" -qb^AC{Ac)"'^ +qb'^C{Ac)"-^ - {^y = 0, l<q<m~l. (17) 

2.3 The case s = 2, m = 3 

The general ideas of the proof can be illustrated for the case where a two stage Runge.-Kutta 
method is applied to problems with a cubic Hamiltonian. Then the abscissae, ci,C2 must be 
those of a third order quadrature rule. Thus, they are the zeros of P2{x) — CPi{x) for some 




Figure 4: This tree has one single leaf to the left and q branches of length two to the right 

C G R, or equivalently the abscissae satisfy the condition 3(ci + C2) — 6ciC2 = 2. We represent 
any 2x2 matrix in in the form 

A = ai_i 15^ + a2a cb^ + ai,2 W^C + ^2^2 cb^C 
There is just one double bush condition p^ with p = l,q = 2, 

b^Ac^ - 2b^CAc = 0. 
Substituting our form of A, and using the quadrature conditions with k < 3, we get 

an + "2,1 + (5 - K) "i'2 + (n - nC) "2,2 = 

where we have made use of the fact that b^c^ ~ i + M^- '^^^ kernel is three dimensional, a 
basis is given as 

A^i = (l-c)6^, N2^21b^ + 3{l~2c)b^C, N3 ^ 2Clb'^ + 3{71 ~ 6c)b^ C 

So any Butcher matrix candidate must be of the form A = cb^ + jiN where iV = wi A^i + ^2^2 + 
V3N3 and the row sum condition (0]) then implies A^l = 0. We obtain after some calculations 
that TV S ker M satisfying this condition must be a multiple of 

iV= ((C-l)l-2Cc)6'^(/-2C) (18) 

There are several possible nonlinear conditions to choose from in order to prove that any candi- 
date solution of the form A = dP' + [3N would require /? = 0. By taking P{x) = G{x) = x{x — 1) 
and R{x) = 1 in p^ and inserting our candidate solution, we find that ^/3^C'^ — which shows 
that one must have /3 = and A = cb^ unless C = 0. 

The remaining case C = can be resolved by using the condition (|17p where upon inserting 
the expression A = dF + jiN one obtains the condition — j^/3^(l + C,Y — Oi therefore /? = also 
for C = 0. 

3 The case of even degree Hamiltonians 

In this section, we prove Theorem [T] for the case that the polynomial Hamiltonian is of even 
degree m = 2s, such that the underlying quadrature is the Gauss-Legendre formula. We use the 
following notation for the standard L^ inner product between functions u and v 

{u,v) — / u{x)v{x)dx 
Jo 



For every non-negative integer q, we let Pq be the Legendre polynomial of degree q 
relative to the interval [0, 1], scaled such that Pq(l) = 1 for every q, and consequently 

(n,p,>^5^ (20) 



(21) 



The polynomials 

Gq{x)^ I Pq-l{t)dt, q>l, 

Jo 
have ioT q > 2 the abscissae of the Gauss-Lobatto quadrature as zeros, and 

The following biorthogonality relations will be useful 

(Gi,P;)-lV^GN, (G,+i,P;)---^V9>2,££N. (22) 

2(?+ 1 

For any quadrature formula, we define the discrete counterpart to the inner product above 

s 

(m, v)d = ^ biu{ci)v{ci) (23) 

and by a slight abuse of language we shall call it the discrete inner product. If the quadrature 
formulas has order m and P and Q are polynomials such that deg P + deg Q < m — I then 

(p,g)z3-(p,g) (24) 

The discrete inner product can still be computed even in cases where it differs from the continuous 
one, the following result which will be of subsequent use, facilitates this in the case where 
deg P + deg Q — m. 

Lemma 7. Suppose that a quadrature rule with abscissae (ci, . . . ,Cs) is exact for all polynomials 
of degree at most m — 1, where s < m < 2s, and let ps = rii=i(^ ~ ^i)- ^^^ ^m ^^ ^ monic 
polynomial of degree m. Then 

(tTto, 1)d = (tT^, 1) - {ps,0rn-s) (25) 

for any monic polynomial 9m-s of degree m — s. 

Proof. Let 6m-i = T^m ^ PsOm-s G ^m-1 for an arbitrary monic polynomial Om-s G Wm-s- Then 

{'^mA)D = {5m-l + Psdm-sA)D = {5m-l, 1)_D = {5m-l, 1) = {T^mA) ~ {Ps,9m~s)- (26) 

D 



We can apply this lemma to obtain the following discrete inner products when Gauss-Legendre 
quadrature is used 



{P2s-r, Pr)D — — 



(G2s-r+l, Pr)D 



l2s-rlr 

7,?(2s + l)' 



li 



£!2 



2s - r + 1 



{P2s-r, Pr)l 



(27) 
(28) 



Here 7^ is the leading coefficient of Pi. 

In analyzing the rank of the linear operator Af , it is useful to work with the transformed 
double bush conditions given in Lemma HI equation p^ . Generally, one may select any suitable 
set of polynomials so that the rank of M is not reduced. In this section we shall make the 
choices Gp and Gq for P and Q, where l<p<q<m — 1. It will also be convenient to write 
the elements A in terms of a basis as follows 






A = 7 , 7 ,ctk,£Ak,e, 
k=l i=\ 



Ak,e ^ Pk-i{c)b'^ PiiC) 



(29) 



The resulting equations for the coefficients ak,i when considering M{A) = are 



^^afe,,((^P-i,^fc-i)D(G„P;)D-(-Pg-i,/'fc-i)i3(Gp,F;)i5)=0, l<p<q<m-l. (30) 
fe=i 1=1 

We prove the following result. 

Lemma 8. Suppose m — Is, s > 2 and c, h are the abscissae and weights of the Gauss-Legendre 
quadrature. Then rank(Af) = 5^ — 1 and ker A/ = span{Ni}. 

Proof. By Lemma [5l clearly rank(Af ) < s^ — 1. Note that the kernel element A^i can be written 
in the format ([2^ as A^i = (1 — c)b^ — -1(^1,1 — ^1,2)- Therefore, it must be true that 
rank(Af) > s^ — 1 if some subset of the conditions (j30l) together with ai.i — cause the 
remaining ak,i to vanish. It is enough to consider just s^ — 1 (linearly independent) conditions 
among the (2s — l)(s — 1). We select the conditions corresponding to 1 < p < g < 2s — 1, and 
such that p + q < 2s + 1. The case s = 4 is reported in figures [5] and H) The {p, q)-e\einent of the 
matrix in Figure [S] corresponds to the condition in (j30p . The numbers refer to the ordering in 
which the conditions are used in the proof. The ones marked {na,nb) are used simultaneously. 
The corresponding ordering of the unknowns a^^i is reported in Figure |6l 
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Figure 5: Ordering of the conditions (p.q). 
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We begin by applying condition (p, s + 1) to Ake for 1 < p < s. Since (P^ applies for all 
inner products, we get from ([20)1 and (|22l) that 



{Pp-l,Pk-l)D {PiGs+l)D - {Ps,Pk-l)D {Pi,Gp)D 



(s + l)(2s + l)(2/c + l) 



Spk^si 



Thus aps = for 1 < p < s. We shall proceed by induction. Suppose it is true that 



aki = 0, 



fc > i + 2, £> i + 1, 



(31) 



which is established for i == s — 2. We now prove, by using conditions {p, i + 2) and {p, 2s — i) 
together with ([5T]) that a^.i = for fc > i + 1, i > i. The first of these conditions applied to Aki 
yields 

{Pp-l, Pk-l)D{Pi,Gi+2)D — {Pi+1, Pk-l)D{Pi,Gp)D 

(HH) applies for all inner products and we conclude, using (PO)) and (1^^ . that condition {p, i + 2) 
implies, after multiplying both sides by {2p — l)(2i + 3) 



- ctp^i+i + ai+2,p-i = 0, p>l 
and for p = 1, multiplying each side by 2i + 3 



(32) 



ai 



,i+i - ^ai+2/ = 0, P=l- 



(33) 



We next consider condition (p, 2s — i) applied to Aki to get 

{Pp-l,Pk-l)D{Pe,G2s-r)D - {P2s-^-uPk-l)D{PiGp)D 

We readily compute {Pp-i,Pk-i)D =- ^ and {PiGp)D = ^^ if P > 1, and (P;,Gi)i3 - 1. 
For (P/, G2s-i)D, (l24t applies when £ < i causing it to vanish by (1221) . For i — i + 1 the total 
degree equals 2s, and by (|28)) 



(^i+l'G'2s-i)£) 



Z+1 

(2s-i) 



TT(^2s-i-l, Pi+l)l 



Nonzero entries for ^ > i + l can be ignored due to the induction hypothesis. Finally, for /c < i + 1, 
one has {P2s-i~i,Pk-i)D = (^2s-i-i,-Pfc-i) = 0. For fc = i + 2, we invoke ((?7)) just to assert 
that (P2S-J-1, Pi+i) D ^ so that this factor can be cancelled in the condition (p, 2s — i) and we 

get 

i + 1 

:OLp,i+i + ai+2,p-i =0, p > 1, (34) 



and 



2s -i 
i + 1 



2s -i 



ai,j:+i 



E 



ai+2,j 



0, p=l. 



(35) 



Combining (j32p and (|M)) we get for p > 1 the system 



-1 1 

■t+i 1 

2s-j -^ 



Q^i+2,p-l 



= 0. 



and thus 



ap,i+i = ai+2,p-i =0, p = 2, . . . , i + 1. 
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The remaining indices to be dealt with in the induction step are ai^^+i and ai+2,i+i and for 
these we consider ([55]) and ([55)) . We use that the only element in the {i + 2)th row {ai+2,e) not 
yet found to be zero is ai+2,i+i and so we obtain also in the case p — 1 a nonsingular 2x2 
system where the two unknowns must satisfy ai^i+i = ai+2,i+i = 0. The induction step is 
completed. The induction proof ensures that all a^ = except possibly ai_i and a2,i, but the 
former is zero by assumption. But the remaining unused condition (p, q) = (1, 2) applied to ^2,1 
yields —1/3 and thus also a2,i = 0. In summary, we have proved that for an s x s-matrix A is 
expressed in the form (P^ with qi^i =0, the conditions ([5(I| imply A — which is equivalent to 
rank(M) > s^ — 1. Combined with the known null- vector A^i = (1 — c)&^ this proves that the 
rank of A is precisely s^ — 1. 

D 

Proof. Theorem [T] (even case). From Lemma HI ([T^ and Lemma [5] we know that any solution to 
the double bush conditions for m = 2s > 3 must be of the form A = cb^ + /3{1 — cjb"^ . But then 
the condition Q immediately implies that 

^1 = c6^1 + /3(l-c)6^1 ^ ;3(l-c) = 

so that /3 = and we are left with the AVF method. D 

4 The case with odd degree 

Suppose now that the degree m of the Hamiltonian is odd. One still needs c and b which satisfy 
the quadrature conditions to order m — 1. This means that it is necessary for the Runge-Kutta 
method to have at least s = (to-|- l)/2 stages such that m < 2s — 1. On the other hand, choosing 
c and b to be such quadrature points and letting A = dP" we have an energy preserving scheme 
with the minimal number of stages. We need to answer whether it is unique. 

The strategy will be the same as in the even case. Now we assume that the quadrature rule 
consists of abscissae (ci, . . . , Cs) which are the zeros of the polynomial Ps — CPs-i for some real ( 
and that bi, . . . ,bs satisfy the quadrature conditions ([3|) for A; < s. The cases C = — 1 and C = 1 
correspond to the Radau I and II formulas in which one has ci = and Cg = 1 respectively. 
Of course C = yields the Gauss-Legendre formula. The case s — 2,m = 3 was discussed in 
Subsection 12.31 and in this section we sometimes assume tacitly that s > 3. 

We let Ri{x) I = 1,2, ... he the polynomials of degree / defined by 

r Pi{x), i^o,...,s-i, 

Ri{x)^ I Ps{x)~CPs-i{x), l = s, 

[ Rs{x)Pi^,{x), l>s + l. 

We consider also the polynomials Fq, defined for every positive integer as 

Fq{x)^ f Rq-l{t)dt, 

Jo 
so Fq — Gq ioi q — 1, ... , s. We observe that for r > s we have 

{P,Rr)D=0. (37) 

for any polynomial P of any degree. 

We will use the following explicit expressions for discrete inner products of polynomials. 
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(36) 



Lemma 9. For any real C we have 

{Pl_,,Fs+r)D =Cj:: f' ^ , Ir-lls-r ^ ^^^^^^ ^ q^ (gg) 

(2s-l)(s + r) 7s_i 

'^•■^•«>°'-^ (2,-l)(. + l) ( (2.-l)|, + 2) ^'-')- <""' 

Proof. The proof relies on Lemma [7] and the three term recursion formulae for Legendre polyno- 
mials, n 

Lemma 10. For all C, ^ —1 there exists a polynomial Pg of degree less than or equal to s such 
that 

(P;,F,+,)z3=0, r = l,...,s-2, (P;,Gi)-0, (41) 

and 

(P;,G2)^0. (42) 

Proof. Consider the vectors fs+r '■= Fs+r{c) G R'*, r = 1, . . . , s — 2 and gi := Gi(c). Then 
dim(span{/5+i,...,/2s-2,5i}) = 7 < s - 1, dim (spanj/^+i, . . . , /2s-2,gi}^") = « - 7 > 1. 

The superscript "J-^)" denotes the complementary vector space with respect to (•, ■)d interpreted 
as an inner product on R'*. Then there exists 5 7^ and g £ spanj/^+i, . . . , /2S-2, .9i}^" • We 
write g by using the basis P{{c), . . . , Ps{c) of R^, i.e. 



Y^veP^ic). 



«=i 



We define Ps := X]f=i^^^i- ^Y construction this polynomial satisfies the orthogonality condi- 
tions of dH]). The condition ^^ is 



(G2,P.)d=^^'^(G2,p;)i, = -^«i, 



and is nonzero if and only if vi ^ 0. We now prove that if C 7^ ^1 then v\ ^ 0. Consider the 
(s — 1) X s matrix V with entries 

r,,, := {F,+,,P'^D, ^ = 1, . . . , s - 2, j = 1, . . . , s, r,_i., := {Gi,P'^d - 1, j - 1, . . . , s. 

We observe that 

r,j=0, i+j<s~l. 

The conditions (|41[) for the vector w := [wi, . . . , Vg] can be written as 

Tw^O. 

Let us define F to be the (s — 1) x (s — 1) Hessenberg matrix whose columns are the last s — 1 
columns of F, and denote by v the vector v := [v2, . ■ ■ ,Vs]'^ , then we have 

Fw = —vi Fei, 
13 



where ei is the first canonical vector in R^. We also note that Tei — ei. If C = 0, f is upper 
triangular because the entries Ti^s-i = {Fs+i,P's_i)D = 0, z = 1, . . . , s — 2, and F is invertible 
because (F^+i, P^'_j_|_;^}£) ^ for i = 1, . . . , s — 2, see Lemmaini 

If <^ 7^ 0, due to the Hessenberg form of F and the fact that Ti^s-i = {Fs+i,P's_i)D ^ 0, one 
concludes that f is invertible if and only if the two last columns of f are linearly independent. 
By Lemma El if C 7^ ^ 1 1 the two determinants 



det 



{Fs+1,Ps-i)d {Fs+i,Ps)d 

1 1 



det 



{F,+2,PLi)d {Fs+2.P's)d 
1 1 



cannot be simultaneously zero. Thus, detF ^ 0, and we can write 

V = —viT~^ ei. 

As a consequence v = [wi, w]^ is a non trivial solution of Fw = giving g = X^fci "^iPgic) j^ 0, if 
and only if wi ^ 0. D 

In the sequel we will use Ps as characterized in Lemma |10l whenever (^ 7^ — 1. When <^ = 1, 
one may take Pg = Ps + Ps-i — 2Pi, and when ( — 0, Pg — P2 — Pi- When C, = —1, wc will use 
instead P^ := Ps - Ps-i, for which only gT]) hold, but not (gg]) as (P^, G2) = 0. 

We express any A in a form similar to (1^^ . 



A= 2^ ak,ePk-i{c)b^Pi{C), 
k,e<s 



(43) 



where we use the notation Pi :— Pi for /<s — 1, ifCT^O, and Pg := Pi for / 7^ 2 and P2 ■= Ps, 
if C = 0. 

Lemma 11. Suppose m — 2s ~ 1, s > 3 and c,b are the abscissae and weights of the Radau 
quadrature with C £ R \ {-!}. Then rank(A'f ) > s^ - 3. // C = -1- then rank(Af ) > s^ - s - 1. 

Proof. For each 1 <p<q<2s — 2 we consider the condition obtained from Lemma [5] by 
choosing P{x) — Fp{x) and Q{x) = Fq{x) as defined in ([55)) . Condition (p, g) applied to A is 



•■p,q 



{A) 



k,e 



Ctkl 



{{Rp-l, Pk-l) D {Pe,Fq)D — {Rq^l,Pk-l)D {PiiFp)d] ■ 



(44) 



Due to the orthogonality properties of the polynomials Rp and Fq, all conditions vanish identically 
iov s + 1 <p < q. We show that if ai^s — 024 = ct2,s = and C 6 R- \ {0, —1}, then mp^q{A) = 
for 1 <p<q<2s — 2 implies that ak/ = for all fc,^ = 1, . . . , s. Analogously, if C = —1, 
ai,s — a2,i = a2.s = and in addition ai_s ~ for £ = 3, . . . , s, then •mp^q{A) = for 
l<p<q<2s — 2 implies ak,i = for all fc, ^ = 1, . . . , s. 

It is enough to consider a subset of linearly independent conditions: s^ — 3, V^ £ R.\ {0, —1}; 
and s^ — s — 1 for C, — —1 respectively. 

These are all the conditions corresponding to {p, q) — (j, s+r), r — 1, . . . , s — 2 and j — 1^ . . . , s 
and successively (p, q) — {j,s — r), r = 0, . . . ,s — j — 1, and j = 1, 2. In figures [7] and H] we show 
the ordering of the conditions and of the unknowns in the case s = 5. The numbers refer to the 
ordering in which the conditions are used in the proof, and in which the ak,e will be shown to 
vanish. 

We start considering the case C 7^ 0, the proof is similar in the case C = and we will highlight 



later the difference. We first consider conditions {p, q) = {j, s + r) and the unknowns a 



j,s-r, 



for 
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Figure 7: Ordering of the conditions (p.q). 
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Figure 8: Ordering of ak,i 



a fixed j and r — 1, . . . , s — 2, and j — 1, . . . ,s. For these conditions, due to ([571) . PH) simplifies 
into 

J2(^3/{PiFs+r)D^0, j^l,...,S. (45) 

/ 

The last term of the sum at the left hand side vanishes due to {Pg,Fs+r)D = 0, see ((4T|) . For 
£ < s — {r + 1) we obtain that 

{P',,Fs+r)D - {PiFs+r) = / P',{t)F,+r[t)dt = 0, 

JO 

this follows using integration by parts and the orthogonality properties of the Legendre polyno- 
mials. So for r = 1 we get the equation 

a,,,_i(F,'_i,F,+i)z,=0, 

implying that aj,s-i = by (j38p . Proceeding by induction over r we similarly obtain that 

aj^s-r{P's-r^ Fs+r)D =0, r = 2, . . . , s - 2, 

and by ([35]), this in turn implies that aj,s~r = respectively for r = 2, . . . , s — 2. 

For ttj^i and j = 3, . . . , s, we consider the conditions p = 1, q = j leading to the equations 

a,u(P{,Gi)=0, 

implying aj^i = 0, since {P{,Gi) = 1. Here we have used that ai^i = for I — 2,...,s, 
{Pi,G,) = 0, for j = 3, . . . , s and (P,', d) = 0. 

For ai^i we use (p, q) — (1, 2), using au = for £ = 2, . . . , s and q;2,£ — for £ = 1, . . . , s we 
simplify the equation into 

ai,i{Pi,G2)^0, 

giving ai^i ~ since (P{, 6*2) 7^ 0. The proof is here concluded for the case C = — 1- 

We next consider conditions {p,q) — (2, s — r), r = 0, . . . , s — 3, leading to the equations 

o X! "2,<?(-Pf , Gs-r)D - , _ s ^ as-r,e{Pe,G2)D = 0, 

since we have shown that a2,£ — for £ = 2, . . . , s — 1 and as-r,i = 0, ^ = 1, . . . , s — 1, and by 
hypothesis, a2,i = «2,s = 0. We are left with the equation 

as-r^s{Ps,G2)=Q. 
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By dHI) (P,', G2) ^ for C 7^ -1, so we conclude a^_r,s = 0, r = 0, . . . , s - 3. 

This concludes the proof, for the cases C G R \ {0}. In the case ( = 0, P{{c), . . . , P^{c) is no 
longer a basis of R* because Pg ^ P2 — Pi ■ We then consider P2 := Ps in P5)) and proceed as 
in the case, CgR\{0,— 1}. The equation (HSI) becomes in this case 

s-l 

"^,2(^2, ^s+r)D + Yl "^jAPi^ Fs+r)D = 0, 
£=s-r+l 

leading to aj^2 ~ for r = 1 and, by induction, to aj^s-r = for 1 < r < s — 2, by using that 
{Ps-r+iJ Ps+r)D 7^ 0. The rcst of the proof is the same as in the previous case. D 

The following two discrete inner products are easy consequences of Lemma [7] 

{Ps+r-l,Ps-r)D = 7^ 7 (46 

{Gs+r,Ps-r)D - '^ {Ps+r-l, Ps-r) D (47) 

s + r 

To establish that the bounds for the rank of M are sharp, we shall simply derive a suitable set 
of linearly independent matrices in the kernel of M, and we make the ansatz that these kernel 
elements are all rank one matrices. We use the generic representation 

s s 

U{c)b'^V{C), U{x)=Y,Uk,Pk-i{x), F(^)=^t.,P;(a;) (48) 

fc=i 1=1 

s 

and for convenience, we shall define vi^ :— ~ \^ vi. 

i.=\ 

Lemma 12. Let s > 2, m = 2s — 1 and suppose that the quadrature rule (c, b) has order 2s — 1 
where Ci, . . . , Cs are the zeros of the polynomial Ps{x) — (^Ps-i{x) for some real (^. Then, if C, i= —1 
there exists a basis for ker M consisting of matrices 

N, = W{c)b^V\C), i^ 1,2,3. 

for polynomials U^ £lls, V^ € Us that can be written in the form 



W{x) = ^4:)p,_i(x), V\x) = ^«f P;(:r) (49) 



fc=l £=1 

having the propertied 



i = l: 


n^^^ = I, 


u^' = -1, 


v'^' = 1, 


i = 2: 


uf' = 1, 


4^) = o, 


^' - 0, 


i = 3 : 


uf' = 0, 


v^ - 1, 


.(^):=^.f) 



= 0. 



If C y^ 0, then one may choose Wg =1, i = 2, 3. 

^The first basis element (i = 1) is nothing else than ANi = 4(1 — c)b^ so that all the coefficients uj^ ,v^ not 
listed are zero. 
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Proof. We first intend to show that the equations 

s 

Y, UkVe{{Pp^uPk-i)D{Gq,P^)D - {Pq-i,Pk-i)D{Gp,P^)D) =0, l<p<q<2s^2, (50) 
k,e=i 

have a solution with the given preset values from the lemma. We already know from Lemma [5] 
that U^{c)b^V^{C) is a kernel element so henceforth we consider only i = 2, 3. 
We treat the case C = separately. Consider 

u(2) = (i,0,-l,0,...,0)^, u(3) = (0,1,-1,..., Of, (51) 

v^'^ = (0,1,0,..., Of, «(3) = (1,-1,...,0)^. (52) 

In these cases, all discrete inner products in (|5(7)) equal the continuous ones for all {p, q) because 
C = implies that the quadrature rule has order 2s, and all conditions are readily verified. 

Until the end of this proof, we assume that C 7^ 0. We next use the conditions {p, q) of ([SO)) , 
setting p = 1,2 and 3 < q < s. We apply ^U^, ([^ and for q = s we use also (|i71) with r = 0. 
We get 

Vp-iUk = UpVk-i, 3 < fc < s - 1, (53) 

Vp-iUs = Up{vs-i - C,Vs), (54) 

The conditions with p = 1 are useful for kernel elements of type i — 2, and those with p — 2 
can be used for the i — i type. We now select the conditions [p, s + r) in ([50)1 where p = 1,2 
and 1 < r < s — 2. In this case, because of ([^0)1 and ([M)l . the discrete inner products will vanish 
whenever i < s — r — 1 and k — l<s — r — 1. We get 

s s 

Up ^ Vt{Gs+r,P'f,)D+Vp-i ^ Uk{Ps+r-l,Pk-l)D =0, p= 1,2. (55) 

£^s — r A;— s* — r+1 

Substitute ([55)1 and (fM|) into these conditions and change summation index 

(56) 

for 1 < r < s — 2 and p ~ 1,2. We have thus derived s — 2 conditions for the s — 1 unknowns 
V2, . . . ,Va. Choosing v^ = 1 there is an upper triangular system to be solved for V2, . ■ . , v^-i. This 
system is nonsingular, since the pivot elements, which can be computed explicitly from (|46p and 
(|Tf)) . are nonzero whenever C 7^ 0. Note that one must choose p — 1 for the kernel element of type 

(2) (3) 

i — 2 and p = 2 for the i = 3 type to have Up ~ 1. So wc conclude that Vg ~ v^ for 2 < £ < s. 
In order to solve for uj! , i — 2,3 from ([55]) - ([5l)) . we need to make sure that Vq = — ^^ w) 7^ 

(3) (2) 

for i = 2, and v^ ^ for i = 3. In the former case, note that Vg = would lead to 

t;^2) ^0, e=l,...,s-2§^, v^^}^ = C dSl, and thereby -v^^^ = v^^}-^ + vP =(+1 = which 

(3) 
is the exceptional case excluded in the lemma. Similarly, v^ ~ would also imply C = — 1. We 

conclude that all parameters u^ , u^ have been determined for ( 7^ 0. It is left to the reader 
to verify that with these choices, all the remaining conditions {p,q), 3<p<q<2s~2 are 
consequently satisfied. 

We already know from LcmmafTTjthat rank(M ) > s^ — 3 for C 7^ — 1, it therefore just remains 
to check that Ni,i — 1,2,3 form a linearly independent set. Suppose 

3 

i=l 
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Multiplying this matrix from the right by the vector c, we obtain 

aiU\c)~a2vl^^U\c) = 0, 

where we have used that {V'^,Gi)d = —Vq = 0. Thus, since Wq ^ for ( =^ —I and since 
clearly U^{c) and C/^(c) are linearly independent we conclude that ai = a2 = so that also 

as = 0. n 

As a by-product of the proof of the previous lemma, we obtain 

Lemma 13. ker A/ consists of matrices A such that rankA < 2. More precisely, for C ^ —1, 
one has 

U'^{x)^U^{x)-U\x), V'^{x)^V^{x)+vfW^ (57) 

We now consider the kernel of M for the exceptional case C = — 1. 

Lemma 14. For ^ = — 1 rank(M) — s^ — s — 1. The following s + 1 matrices constitute a basis 
for ker M 

iVi=4(l-c)6^, iV,+i=P,_i(c)6^(P;(C)-P;_i(C)), i = l,...,s. 

Proof. The linear independence of these matrices is evident. The identity x{Pl — ^s'_i) = s{Ps + 
Ps~i) shows that {Gq,Pg — P's^i)d = for all q > 1 so that ([50]) is trivially satisfied for all p 
and q. Since there are s + 1 linarly independent elements, Lemma 1 1 1 1 ensures the rank of M is 
precisely s^ — s — 1 and that the given set of matrices forms a basis for ker M when C, — —1. D 

At this point we know that a necessary condition for a Runge-Kutta method to be energy 
preserving for polynomial Hamiltonians is that its Butcher matrix is of the form A — cb^ + N 
where N £ kerM. Now, we make use of the row sum condition ^ to infer the additional 
requirement on N that N ■ 1 = 0. The following lemma is easily proved. 

Lemma 15. For all (^ G R, the intersection o/ker Af and the set of all s x s matrices with zero 
row sum, i.e. N -1 — 0, is a one dimensional subspace ofR^^^ consisting of matrices of rank at 
most one. 

From the preceding lemma, we now conclude that any Butcher matrix we search for is of 
the form cb^ + (3 N where N — Ub^V G ker Af. We show that such candidates are indeed 
incompatible with the nonlinear triple bush conditions presented in subsection 12 . 2 1 unless /? = 0. 
Of particular use to us is the condition obtained by choosing P{x) = Q{x) — Gp{x) and R{x) = 
Gi{x) = a; in (|16l) . Because of the symmetry, and because Gp{l) = Spi, we get the following 
simple special case: 

2b'^Pp^i{G)ACAGp{c) ~ 2Spib'^CAGp{c) - b'^{AGp{c) AGp{c)) + i^l=0, (58) 

where i^i = 2, ^^2 = g and Vp — when p > 2. 

Lemma 16. If A = cb^ + /S N is a solution to the triple bush condition (|58p where N — 
U{c)b^V[C) e ker Af with ( ^ {-1, 0}, then ^ = 0. 



Proof. We substitute such a solution into (ISSp with p = 1, and p = 2. The terms coming from 
cb'^ cancel since the AVF method is energy preserving. The linear terms in (3 are proportional to 
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U2Vo — uivi and therefore vanish by (j50p using p = 1, q = 2. The quadratic terms are computed 
by means of (|^D|) and (|22p . and we get the condition 

{2upxV + Vp^iU,Vp-iU)D — 0, P = li2. 



(2p-l)2 
We invoke ^3^, ([Ml) to substitute for Vp^iU 

fl2 






(2.tF + V^,^)z5=0, where V^} Vk-iPk-i - (vsPs-i 



The inner product can be worked out, and one finahy gets the condition 

^ ■w2«2(C+ 1)2 = 0, p=l,2. 



(2p-l)2 P '' 

From Lemma \n\ and Lemma [T^] we deduce that ui = and U2 = simuhaneously is impossible 
for any rank one kernel element. For all C 7^ 0, u^ = 1, so the lemma holds as stated. D 

We now consider the case C, — ^ which corresponds to the Gauss-Legendre quadrature for- 
mula. From (|52l) we see that the second kernel element has V'^{x) — P2{x) which means that 
U^{c)b'^V^{C)l — U^{c){P2,1)d — 0, this causes the one dimensional subspace of kernel ele- 
ments satisfying the row sum condition to be the span of 

N - U{c)b'^V{C) = (Po(c) - P2{c))b'^P^{C) (59) 

where we have skipped the superscripts on U and V for ease of notation. We use the condition 
((T7)) with q — s which reads, after inserting A — dF + fiN , and using ([/, 1)d ~ (V, x)£i = 1, 

{{^x + puy, i)d + § (x, {\x + iiuy~^)D - P{xv, {\x + puy'^)D - {\y = o. (60) 

We observe that these inner products involve polynomials of degree 2s whereas the underlying 
quadrature formula for this case has order 2s and is therefore exact for all polynomials of degree 
at most 2s — 1. Clearly, the condition (j60l) would hold exactly for all /3 if all discrete inner 
products were replaced by continuous ones. We therefore conclude that it is only necessary to 
retain terms arising from the leading order 2s, these are the terms multiplying /?*. This results 
in the condition 

Is 

and we have proved that any solution of the form A = dP^ + f3N requires /3 — 0. 

The final exceptional case is C = ~1 which corresponds to the Radau I quadrature and ci = 0. 
The one dimensional subspace of elements in ker M which satisfy the row sum condition is in 
this case given as the span of the matrix 

N = U{c)b^V{C), U{x) = Po(a;) - Pi{x), V{x) = {-lYP[{x) - P'^_^{x) + P^{x) 

We substitute A = cb^ + /3N into the triple bush condition (|T6| with P{x) = Q{x) = xG2{x), 
R{x) = 1. After evaluating all inner products, we find that the condition yields 

- l^' = 0. (61) 

and therefore we must have /3 = and the only possible energy preserving method with ^ = — 1 
is the AVF method, A — cb^ . We summarize these findings 
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Lemma 17. Let ( G { — 1,0}, N e kerM and consider the matrix A = dF + fiN , where 
N ■ 1 = 0. For the Runge-Kutta method with coefficients {c,h,A) to he energy preserving, one 
must have j3 — Q. 

Proof. Theorem [T] (odd case) . By Lemma [16] and Lem.ma [17] the theorem, is proved. 

D 

Acknowledgments. 

The first two authors would like to acknowledge the support from the IRSES project CRISP, and 
part of the work was carried out while the authors were visiting Massey University, Palmerston 
North, New Zealand and Latrobe University, Melbourne, Australia. 

References 

H. Berland and B. Owren. Algebraic structures on ordered rooted trees and their significance 
to Lie group integrators. In Group theory and numerical analysis, volume 39 of CRM Proc. 
Lecture Notes, pages 49-63. Amer. Math. Soc, Providence, RI, 2005. 

J. C. Butcher. Numerical methods for ordinary differential equations. John Wiley & Sons 
Ltd., Chichester, second edition, 2008. 

E. Celledoni, R. I. McLachlan, D. I. McLaren, B. Owren, G.R.W. Quispel, and W. Wright. 
Energy-preserving Runge-Kutta methods. ESAIM: M2AN, 43(4):645-650, 2009. 

E. Celledoni, R. I. McLachlan, B. Owren, and G.R.W. Quispel. Energy-preserving integra- 
tors and the structure of B-seies. J. of FoCM, 10:673-693, 2010. 

P. Chartier, E. Faou, and A. Murua. An algebraic approach to invariant preserving inte- 
grators: the case of quadratic and Hamiltonian invariants. Numer. Math., 103(4):575-590, 
2006. 

E. Faou, E. Hairer, and T. L. Pham. Energy conservation with non-symplectic methods: 
examples and counter-examples. BIT, 44(4):699-709, 2004. 

E. Hairer. Energy-preserving variant of collocation methods. Journal of Numerical Analysis, 
Industrial and Applied Mathematics, 5(l-2):73-84, 2010. 



[s; 



[lo; 



[11 



E. Hairer, C. Lubich, and G. Wanner. Geometric numerical integration, volume 31 of 
Springer Series in Gomputational Mathematics. Springer- Verlag, Berlin, second edition, 
2006. Structure-preserving algorithms for ordinary differential equations. 

E. Hairer, S. P. N0rsett, and G. Wanner. Solving ordinary differential equations. I, volume 8 
of Springer Series in Computational Mathematics. Springer- Verlag, Berlin, second edition, 
1993. Nonstiff problems. 

F. lavernaro and D. Trigiante. High-order symmetric schemes for the energy conservation 
of polynomial Hamiltonian problems. JNAIAM, 4(l-2):87-101, 2009. 

B. Leimkuhler and S. Reich. Simulating Hamiltonian dynamics, volume 14 of Gambridge 
Monographs on Applied and Computational Mathematics. Cambridge University Press, 
Cambridge, 2004. 

20 



[12] B. Owren and A. Marthinscn. Runge-Kutta methods adapted to manifolds and based on 
rigid frames. BIT, 39(1):116-142, 1999. 

[13] G. R. W. Quispel and D. I. McLaren. A new class of energy-preserving numerical integration 
methods. J. Phys. A, 41(4):045206, 7, 2008. 



21 



